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Abstract 


The 3DHZETRN code, with improved neutron and light ion (Z < 2) transport procedures, was 
recently developed and compared to Monte Carlo (MC) simulations using simplified spherical 
geometries. It was shown that 3DHZETRN agrees with the MC codes to the extent they agree with 
each other. In the present report, the 3DHZETRN code is extended to enable analysis in general 
combinatorial geometry. A more complex shielding structure with internal parts surrounding a 
tissue sphere is considered and compared against MC simulations. It is shown that even in the 
more complex geometry, 3DHZETRN agrees well with the MC codes and maintains a high degree 
of computational efficiency. 

Introduction 

A 3D version of HZETRN has recently been developed and tested in simple geometries and compared to 
available Monte Carlo (MC) codes [Wilson et al. 2014a-c, 2015]. In the first study [Wilson et al. 2014a, b], an 
aluminum sphere with a radial thickness of 20 g/cm 2 was examined, and comparisons were made between 
3DHZETRN results and MC simulations at locations within the sphere exposed to solar particle event (SPE) and 
galactic cosmic ray (GCR) boundary conditions. Well defined convergence tests were performed that spanned the 
transport formalism, covering the straight ahead approximation (N = 1), bi-directional approximation (N = 2), and 
more detailed 3D treatment (N > 2) for neutrons and light ions. Note that in the formalism of Wilson et al. [2014a-c, 
2015], N denotes the number of transport directions used to describe the assumed isotropic neutron and light ion 
fluxes at lower energies. The convergence tests and comparison to MC results established improved description of 
the neutron and light ion solutions within the 3D formalism. It was also shown that although there are still noticeable 
differences associated with the nuclear production models used in the codes, 3DHZETRN agreed with the MC 
models to the extent they agree with each other in finite geometry. The main difference between the transport codes 
was found to be in the computational costs. While 3DHZETRN results were generated in several seconds on a single 
CPU, the MC code run times were several orders of magnitude larger (-10 8 seconds) and required high performance 
computing clusters. 

In the next study [Wilson et al. 2014c, 2015], a tissue sphere with a radial thickness of 15 g/cm 2 , defined by 
the International Commission on Radiation Units and Measurements (ICRU) [ICRU 1993], was surrounded by an 
aluminum spherical shell with thickness 20 g/cm 2 , and similar convergence tests and comparisons to MC were 
provided. Improvements in the low energy neutron and light ion spectra as a result of 3D corrections were again 
clearly established. Despite the added geometric complexity of the two material sphere configuration, agreement 
between all the codes was actually improved compared to the previous single material aluminum sphere. The 
improved agreement was attributed to the hydrogen content of the ICRU tissue sphere. For energies below -100 
MeV, elastic collisions between neutrons and hydrogen dominate the neutron transport processes. These elastic 
collisions, on average, transfer half the neutron energy to the target hydrogen, which in turn, very rapidly attenuates 
the neutron energy spectrum. Although neutron production cross sections show significant variation amongst the 
codes for aluminum targets, as seen in the first study [Wilson et al. 2014a,b], neutron-hydrogen elastic cross sections 
are more precisely and accurately represented in the codes either through evaluated nuclear data files or detailed 
parameterizations [Wilson et al. 1991]. Again, it was found that 3DHZETRN was in good agreement with the MC 
simulation results, and the primary difference between the codes lay in the associated computational costs, where 
3DHZETRN is several orders of magnitude faster. 

Herein, a more complex and realistic shielding geometry with multiple parts is considered. An ICRU tissue 
sphere is surrounded by a cylindrical aluminum shell containing two additional internal aluminum boxes. The sizes 
and locations of the internal parts are chosen to provide partial shielding of the ICRU sphere, thereby emphasizing 
3D effects in the transport process. The chosen geometry is sufficiently complex to further test 3DHZETRN while 
still allowing engagement of various MC codes to verify the solution methods. 

In this report, transport code development efforts will be briefly reviewed with an emphasis on the most 
recent extensions leading to a generalized version of 3DHZETRN applicable in combinatorial geometry. As in the 
previous studies [Wilson et al. 2014a-c, 2015], convergence tests are performed and comparison to MC simulations 
are provided. The next step in 3DHZETRN development will allow the use of more complex and realistic geometric 
models, including human phantoms, so that simple mapping of the present methodology into more realistic 
applications can be studied. Although a final solution to engineering design problems is not yet at hand, the current 
status of deterministic methods agrees with MC codes to the degree that various MC codes agree among themselves 
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in many cases. The main limitation within 3DHZETRN remains in the nuclear databases and requires additional 
experimental measurements and nuclear model development. 

3DHZETRN 

The 3DHZETRN theoretical formalism has been provided in prior reports [Wilson et al. 2014a-c, 2015] 
and will not be repeated in detail here. However, an overview is given here to provide clarity in the notation and 
terminology used later in this paper. Additional discussion is given in this section regarding extensions of the 
present code to general combinatorial geometry. 

Theoretical Formalism Overview 

The linear Boltzmann transport equation within the continuous slowing down approximation for the flux 
(or fluence) density, 0 .(aj, ft, E ) , of a j type particle is given by [Wilson et al. 1991, 2005] 


b [4>j(x,n,E)] = J2fzf 4 ^jk( E ’ E '’n,n')<fr k ( x W, E ')dn' dE ', (i) 

k 

where the differential operator on the left hand side is defined as 

B [</>. (x, Cl,E)} = Cl- V(f>j (as, n, E )~Y~^l S j (a, Cl, E ) ] + a. (Ety (x, Cl, E) . (2) 

In equations (1) and (2), Aj is the atomic mass of a type j particle, Sj(E) is the stopping power of a type j ion with 
kinetic energy E (vanishes for neutrons), Oj(E) is the total macroscopic cross section for a type j particle with kinetic 

energy E , and a - k {E,E\ft,ft } ) is the double differential macroscopic production cross section for interactions in 

which a type k particle with kinetic energy E and direction Ft ' produce a type j particle with kinetic energy E and 
direction ft . 

Solution methods are developed by separating the double differential cross section for neutron (j = n ) 
production and the particle fluxes into forward and isotropic components. The forward components are associated 
mainly with higher energy direct quasi-elastic events and projectile fragmentation products [Wilson 1977, Wilson et 
al. 1988], and the isotropic components are associated with lower energy secondary particles, including target 
fragments. 

The forward component is first solved within the straight ahead approximation, wherein all particles are 
assumed to travel along a common axis (12 ~ O'). This allows previously developed and highly efficient numerical 
solution techniques to be utilized [Wilson et al. 1991, Slaba et al. 2010a]. The forward solution then gives rise to a 
source of isotropically produced neutrons, given as 

tn,iso&n,ci 0 ,E) = jy nK UE,E\n,n 0 )<i )kJor {x,ci Q ,E')dE' , O) 

k 

where (j n k,iso(E,E\Fl,Fl') is the isotropic component of the neutron production cross section, and <\ \j or (x,fl,E') is the 
forward component of the flux. The symbol Qo denotes the direction of the inbound forward flux arriving at location 
x, and ft denotes the transport direction within the geometry along which the isotropic flux is being computed, as in 
Fig. 1. The transport equation for the isotropic neutron flux is given by 


[«*V + a n (E)}^Jx,n,E) = f™f^a nn (E,E\n,n')^ i Jx,n\E')dn'dE' 


0 ’^) 


( 4 ) 


2 



12o 



Fig. 1. Geometry depicting relationship between the forward direction 12o and transport direction 12 at location*. 


Equation (4) is solved along the transport direction 12 within the bi-directional approximation discussed in 
detail elsewhere [Slaba et al. 2010b]. A final step in the solution methodology is to compute the source of light ions 
produced from the lower energy isotropic neutrons. Once this source is computed, the isotropic component of the 
light ion flux is solved under the assumption that no further nuclear collisions occur, giving partial 3D treatment to 
low energy charged particles. 

A remaining detail in the transport formalism is to define a discrete number, N, of transport directions, 12, 
over which the isotropic flux is evaluated. In connecting to historical code development, the N = 1 solution 
corresponds to the usual straight ahead approximation wherein all particles are assumed to follow a common axis. In 
this case, there is only one transport direction 12, and 12 = 12o. The N = 2 solution corresponds to the bi-directional 
transport approximation, wherein particles are allowed to propagate either straight forward or straight backward. In 
this case, there are two transport directions (12 1 = 12o and Q 2 = -!2o). Values of N > 2, correspond to a more detailed 
three dimensional (3D) description of the isotropic flux component and are the emphasis of recent advances in code 
development [Wilson et al 2014a-c, 2015]. In those previous studies, values of N = 1, 2, 6, 10, 14, 18, 22 were 
considered, as shown in Fig. 2. 



Fig. 2. Transport directions used to describe isotropic component of flux. 


These distributions were developed by prescribing a number of latitudinal (No) and longitudinal (Np) 
directions along with the two polar directions. For the present study, a generalized algorithm for computing the 
transport direction vectors, 12, was implemented into the code to enable more comprehensive testing in complex 
geometry, as will be shown in the next section. 

The algorithm is described as follows. The number of latitudinal (No) and longitudinal (Np) directions is 
chosen. This gives the total number of directions 


N = N a Np + 2 , (5) 

where the two additional polar directions (forward/backward) have been added. The longitudinal angles (/?) and 
latitudinal angles (at) are computed as 
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i — 1 


( 6 ) 


A = 2?r - 


iV„ 


a- = arccos 


a, + —(2i -1) 


( 7 ) 


where a = — 1 + 2 / N , 6 = 1 — 2 / N , and A = (b — a) / N a . The directional components are then computed 
using the typical spherical coordinate transformations 


a = 


cos(/3 i )sin(a i ) 

sin(/3 i )sin(a i ) 

cos(aJ 


( 8 ) 


In the following section, convergence in the induced radiation field within the ICRU tissue sphere is examined as a 
function of N, and transport direction distributions as large as N= 86 are considered. 

Combinatorial Geometry 

The present work required extensions to 3DHZETRN to enable general combinatorial geometry to be 
studied. As a result, an efficient computational ray-tracer was developed and integrated into the code. 
Fundamentally, the ray-tracer provides the thickness and type of each material traversed in moving from any 
location in space along an arbitrary direction through the geometry. Such information is required in two places 
within the 3DHZETRN formalism (see Fig. 1). First, to evaluate the forward flux along the direction Qo, the ray- 
tracer provides the ordering, material types, and thicknesses traversed in passing from the outermost geometric 
boundary to the point of interest, x. Second, to evaluate the isotropic flux along the transport direction £1 passing 
through x , the ray-tracer is called again to provide similar information. 

Solid spheres, boxes, cylinders, and ellipsoids can be easily ray-traced since the equations describing 
intersection between any ray in space and one of the solids have a closed form analytic solution. Additional 
geometric and computational complexity is introduced when these simple solids are nested inside each other but 
must be included as a capability. 

The currently implemented ray-tracer allows any number of the simple solids to be nested inside each 
other, allowing relatively complex geometries to be defined. Partially nested solids (i.e. solids that are not entirely 
contained by another solid) are not currently supported. Nesting geometric solids is achieved computationally by 
requiring simple containment (mother-daughter) relationships to be explicitly defined prior to execution. In 
particular, the definition of each object prior to execution includes a field identifying its immediate mother volume 
(i.e. the solid that contains it). From this basic information, a complete mother-daughter tree is generated within the 
code for the ray-tracer upon execution. 

The overall procedure of the ray-tracer is described as follows. For a given point and ray-trace direction, 
intersection distances between the point and the boundaries of each solid in the geometry are computed. If no 
intersection with a given solid occurs, the evaluated distance is set as a negative number. Next, the list of 
intersection distances is sorted from smallest to largest. This sets the ordering of objects and materials traversed 
along the ray. Negative numbers in the list, corresponding to solids that were not intersected, are discarded. The 
remaining list of intersection distances uniquely defines the ordering, material types, and lengths traversed by the 
ray and is then passed from the ray-tracer to the transport procedure. 

The present status of the ray-tracer allows relatively complex shielding geometries to be efficiently studied 
within the 3DHZETRN transport code. Future work will describe implementation of detailed human phantom 
models into the code. Also, extensions to the ray-tracer to enable compatibility with shielding thickness 
distributions, as generated in the computer aided design vehicle optimization analyses, are forthcoming. Both of 
these advancements will provide the final connections between the computationally efficient 3DHZETRN code to 
engineering design and radiation risk studies for space applications. 
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Complex Inhomogeneous Geometry 


As in prior studies, geometry of sufficient simplicity will be chosen to allow MC simulation using available 
computer resources in a reasonable commitment of time but of sufficient complexity to allow testing of the 
3DHZETRN procedures. The geometry under consideration is shown in Fig. 3. A 30 cm diameter tissue equivalent 
sphere [ICRU 1993] is placed within a 2.47 cm (5 g/cm 2 ) thick cylindrical aluminum shell (outer diameter 75 cm, 
outer length 175 cm) with 4.9 cm (10 g/cm 2 ) thick end caps. In addition to the ICRU sphere at the center of the 
aluminum cylindrical shell, an interior aluminum box with dimensions 46.3x14.8x14.8 cm (upper right comer in 
Fig. 3) and a second aluminum box with dimensions 22.2x22.2x37.0 cm (lower left corner in Fig. 3) are also placed 
inside the cylinder. The dimensions of these boxes in units of g/cm 2 are 125x40x40 and 60x60x100, respectively. 
Target points (detectors) are placed within the ICRU tissue sphere, similar to the previous study [Wilson et al. 
2014c, 2015]. The detector locations appearing at the top and bottom of the ICRU tissue sphere will be referred to in 
subsequent figures, tables, and discussion as z = 0 and z = 30, respectively. Additional target points within the ICRU 
tissue sphere are therefore defined by their location between the top and bottom detectors. For example, the detector 
location in the center of the sphere is referred to as z = 15. 

The external radiation environment (external source boundary condition) is positioned above the geometry 
and directed uniformly downward along Qo as indicated by the arrows in Fig. 3. For all calculations, the boundary 
condition, or external source, is the Webber representation of the February 23, 1956 solar particle event (SPE) 
[Webber 1966]. The Webber 1956 SPE is given by an exponential momentum spectmm with momentum parameter 
po = 100 MV and 10 9 protons/cm 2 above 30 MeV. 

An important feature of the complex geometry is the partial shielding of the ICRU tissue sphere provided 
by the internal aluminum boxes. In particular, the box in the upper right comer casts a shadow onto part of the tissue 
sphere as shown in Fig. 4. Note that similar effects were also shown in a prior study [Wilson et al. 2014c, 2015] for 
various external boundary conditions but in simpler geometry. Fig. 4 shows the source of neutrons induced by the 
Webber SPE event within the complex geometry computed with 3DHZETRN. The rapid attenuation of the neutron 
source intensity along the direction of the incident SPE protons, Oo, is evident, as is the shadow cast by the upper 
right aluminum box onto the tissue sphere. The added complexity of the non-uniformly shielded and partially 
shadowed tissue sphere should accentuate 3D transport processes and provide a more rigorous test for the 
3DHZETRN model. 

The neutron fluence at a point within the sphere depends not only on the distance to the boundary along a 
given transport direction, £1, but also on the intensity of the source along £1 in reaching the boundary. The most 
intense regions of the neutron source must be adequately represented in the domain of the numerical solution. This 
requires sufficient fidelity in the set of transport directions, Q, as defined by equations (5) - (8) in the previous 
section. The convergence rate for increasing number of N will be tested in the following section. 


1956 Webber SPE 




Fig. 3. Complex geometry used in the current study: ICRU tissue sphere shielded by an aluminum cylindrical shell 

with internal parts. 


5 


1956 Webber SPE 
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Fig. 4. Integral (E> 1 MeV) isotropic neutron source (particles/(g-event)) induced by the Webber SPE incident on 

an ICRU sphere shielded by a complex aluminum structure. 


Preliminary testing 


A convergence study was made to examine the solution dependence on the number of longitudinal nodes, 
Np, and latitudinal nodes, N a , used to define the set of transport directions, Q. For this convergence test, the quantity 
studied is dose equivalent at the bottom of the ICRU sphere, since it is typically a more sensitive measure than dose. 
The solution with N = 86 (Np = 12, AU = 7) is taken as the converged solution, and the relative error over a range of 
values is shown in Fig. 5. 

By comparing recent studies [Wilson et al 2014a-c, 2015], it can be seen that the convergence behavior 
depends on the shape and complexity of the shield geometry, as would be expected. For example, in those prior 
studies, results computed in spherical shielding converged well with N = 18 and above. In the present set of results, 
it can be seen that relative error less than 3% is achieved for N a > 4 and Np> 5. 


A T = 44 


N= 58 


N = 12 


N= 86 


Relative error (%) 



Number of longitudinal nodes, 

Fig. 5. Relative error (%) in dose equivalent at the z = 30 detector location (bottom of the ICRU sphere). 
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Also of interest is the dose and dose equivalent at locations within the tissue sphere. Dosimetric quantities 
are evaluated at depths in the tissue sphere at z = {0, 0.007, 0.3, 1, 5, 25, 29, 29.7, 29.993, 30} g/cm 2 with results in 
Tables 1 and 2 for A = 1, 2, 6, 10, 26, 38, 58, 86. The doses near the top and bottom of the tissue sphere are seen to 
converge rapidly with increasing A. The local fields in these locations are dominated by transport through the 
aluminum shield giving improved convergence behavior over the case concerning the surface of a homogeneous 
aluminum shield [Wilson et al. 20 14a, b]. The rate of convergence of dose equivalent is somewhat slower, since the 
details on the penetration of charged particles produced by neutrons are more dependent on angular factors. Still, the 
rate of convergence is quite rapid in the complex geometry with adequate results in only 20 or more transport 
directions (see Fig. 5), similar to the spherical shell shield studied in earlier papers [Wilson et al. 2014c,2015]. 

The solution methodology appears to be nearly converged for A > 20. This is supported in Tables 1 and 2 
by comparing exposure values for different A at a fixed depth. The relative difference, compared to the A = 86 
solution, is found to be less than 1% for dose and less than 2% for dose equivalent for A > 10 at all depths. The 
induced particle fluence is shown in Fig. 6 at the near (z = 0) and distal (z = 30) interfaces. Also shown are the 
lowest order A = 1 and A = 2 solutions, corresponding to the straight ahead and bi-directional approximations, 
respectively. The A = 1 solution has no contributions from the backward leaking neutrons that have been moderated 
within the tissue sphere. This can be seen in the incorrect neutron spectral shape below 20 MeV at z = 0. This error 
is largely corrected by the N= 2 bi-directional approximation. All solutions are nearly converged for A > 10 (i.e. 
within fractions of a percent of the A = 86 result). 

Table 1. Webber SPE dose as a function of number of transport directions, A, at various depths in the tissue sphere 
incl uding those defined by ICRU. 


z (g/cm 2 ) 




Dose (cGy/event) 




A= 1 

N=2 

VO 

II 

A= 10 

A= 26 

A= 38 

ii 

Ua 

OO 

vo 

00 

II 

0 

42.76 

42.78 

42.76 

42.80 

42.79 

42.79 

42.78 

42.78 

0.007 

42.46 

42.48 

42.47 

42.49 

42.49 

42.48 

42.48 

42.48 

0.3 

38.68 

38.76 

38.74 

38.76 

38.75 

38.75 

38.74 

38.75 

1 

32.32 

32.40 

32.38 

32.40 

32.39 

32.38 

32.38 

32.38 

5 

14.56 

14.61 

14.59 

14.61 

14.59 

14.59 

14.59 

14.59 

25 

1.55 

1.55 

1.54 

1.54 

1.55 

1.54 

1.55 

1.54 

29 

1.13 

1.13 

1.13 

1.12 

1.13 

1.13 

1.13 

1.13 

29.7 

1.08 

1.07 

1.08 

1.06 

1.07 

1.07 

1.07 

1.07 

29.993 

1.05 

1.05 

1.05 

1.04 

1.05 

1.05 

1.05 

1.05 

30 

1.05 

1.05 

1.03 

1.03 

1.04 

1.04 

1.04 

1.04 


Table 2. Webber SPE dose equivalent as a function of number of transport directions, A, at various depths in the 
tiss ue sphere including those defined by the ICRU. 


z (g/cm 2 ) 



Dose Equivalent (cSv/event) 



A= 1 

N=2 

A= 6 

A= 10 

A= 26 

A= 38 

S* 

II 

Lh 

OO 

5? 

II 

00 

ov 

0 

65.90 

65.99 

65.69 

66.08 

65.99 

65.94 

65.86 

65.90 

0.007 

62.51 

62.56 

62.37 

62.62 

62.51 

62.44 

62.41 

62.45 

0.3 

54.92 

55.49 

55.26 

55.47 

55.36 

55.28 

55.26 

55.31 

1 

45.73 

46.27 

46.06 

46.24 

46.15 

46.06 

46.04 

46.09 

5 

20.68 

21.03 

20.81 

20.94 

20.81 

20.76 

20.78 

20.82 

25 

2.64 

2.63 

2.56 

2.47 

2.59 

2.57 

2.58 

2.57 

29 

2.03 

2.00 

2.04 

1.90 

2.04 

2.02 

2.05 

2.05 

29.7 

1.95 

1.91 

2.00 

1.83 

1.98 

1.96 

1.99 

1.99 

29.993 

1.91 

1.86 

1.99 

1.79 

1.94 

1.91 

1.96 

1.95 

30 

1.91 

1.84 

1.80 

1.74 

1.85 

1.83 

1.84 

1.84 


The fluence spectra in the neighborhood of the near and distal interface are shown in Fig. 7 for A = 26. A 
rapid transition occurs in passing from aluminum to tissue in the near interface and in passing from tissue to 
aluminum in the distal interface. This can be seen in the proton and 4 He spectra in Fig. 7 below 10 MeV/n. The 4 He 
fluence spectra are shown in the neighborhood of the near and distal interfaces on an expanded scale in Fig. 8. The 
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4 He ions produced in the aluminum shield rapidly form a new equilibrium spectrum within the first few millimeters 
of tissue as observed in the near interface. Similarly, at the distal interface, the leakage of 4 He from the aluminum 
into the tissue shield quickly forms a new, less dynamic, approach to an equilibrium spectrum within the first few 
millimeters of tissue. 

One difference for the present complex inhomogeneous geometry compared to the tissue sphere in close 
contact with the aluminum spherical shell [Wilson et al. 2014c, 2015] is that the aluminum surface in the present 
case is exposed to the secondary light ions over a restricted solid angle whereas the spherical shell configuration is 
over the full hemisphere. This results in some differences in behavior in the transition region. This is seen most 
clearly by comparing the present result in Fig. 8 with the results derived previously for the inhomogeneous sphere 
[Wilson etal. 2014c, 2015]. 

It is found that the details of the 4 He spectra (for example, see Fig. 8) cannot be verified by MC simulations 
since the computational costs required to control statistical fluctuations at low energies are too large. Even the hope 
of defining the transition regions for neutrons and protons has not been successful in MC simulations. Biasing 
techniques could be implemented to resolve these issues in future work. As a result, the larger features generated 
from MC simulations are used to verify the 3DHZETRN results. If the trends in the 3DHZETRN solution over 
larger distances are verified, then it is reasonably assumed that the interface details revealed by 3DHZETRN are 
reasonably correct in the context of the current interaction models. The N = 26 will be used in the MC benchmarks 
discussed in the next section. 




Fig. 6. Webber SPE convergence test at the top (left plot) and bottom (right plot) of a 30 g/cm 2 diameter tissue 

sphere within a complex shield of aluminum. 




Fig. 7. Webber SPE induced fluence within the ICRU sphere at locations near the top (left plot) and bottom (right 

plot) tissue/shield interface with N= 26. 
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Fig. 8. Webber SPE induced 4 He fluence within the ICRU sphere at locations near the top (left plot) and bottom 

(right plot) tissue shield interface with N= 26. 


Monte Carlo Benchmarks 

Using the external source and geometry defined in Fig. 3, the induced neutron and proton fluence spectra 
were evaluated at detector locations within the ICRU sphere by three MC codes (Geant4 [Agostinelli et al. 2003], 
FLUKA [Fasso et al. 2005, Battistoni et al. 2007], and PHITS [Sato et al. 2006, 2013]) and 3DHZETRN. Results are 
shown in Fig 9. The scoring regions (detectors) and physics options used in this set of simulations are identical to 
what was used in Wilson et al. [2014c, 2015] and will not be repeated here. 

The proton fluence predictions among the three codes are in reasonable agreement above 1 MeV with very 
large statistical uncertainty below 1 MeV (not shown). Some of the differences in proton spectra at z = 30 may be 
associated with the simplified assumptions of the straight ahead approximation in 3DHZETRN light ion transport. 
The larger discrepancy lies with the neutron spectra, where disagreement among the MC results brackets the 
3DHZETRN solution. The FLUKA and 3DHZETRN results exhibit a similar inflection in the neutron spectrum near 
10 to 60 MeV that may be related to the use of an early FLUKA approximation to the secondary neutron spectrum 
used in HZETRN [Ranft 1980, Wilson et al. 1988] and carried over into 3DHZETRN. 

To further quantify differences among the codes, the neutron spectra are converted into effective dose using 
conversion factors for isotropic neutrons evaluated by Pelliccioni [2000] using FLUKA. The effective dose on the 
top and bottom of the tissue sphere is given in Table 3. The 3DHZETRN (N= 26) compares well with the three MC 
codes, and is within 5% of the PHITS results. The older HZETRN (N= 1) provides a conservative overestimate as 
expected. An added measure of degree of agreement is to evaluate the relative root meant square (RMS) differences 
in the neutron fluence from the four codes as given in Table 4 on the top and the bottom of the ICRU sphere. Again, 
it is seen that the agreement of 3DHZETRN with the MC codes is about the same as the agreement of the MC codes 
among themselves. 

A principle difference in the three codes is the computational efficiency. The time required by the three codes 
to run these benchmarks is given in CPU seconds used and listed in the first column of Table 5. Although it is 
difficult to favor one code over another on the basis of the veracity of the solutions obtained, there are vast 
differences in computational time required to evaluate fluence spectra even in this simplest of complex geometries 
(tissue sphere within a complex shield of a single elemental material). Still, Monte Carlo methods provide an 
important verification for developing codes capable of meeting operational and especially design requirements. 
Computational speed has proven to be important in spaceflight validation in low Earth orbit (LEO) where the time 
structure of the environment can be used to test various environmental components as was done using Liulin 
instrument measurements onboard the International Space Station (ISS) [Wilson et al. 2007, Slaba et al. 2013]. The 
main limitation on such studies remain the uncertainty in the environmental models (especially for the trapped 
environment) and uncertainty in nuclear cross sections. 
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Fig. 9. Webber SPE benchmark with N= 26 compared with Geant4, FLUKA, and PHITS at detector locations near 
the top (left plot) and bottom (right plot) of the ICRU tissue sphere. 


Table 3. Neutron effective dose values (mSv/event) at z = 0 g/cm 2 (top of ICRU sphere) and z = 30 g/cm 2 (bottom of 
ICRU sphere). Propagated MC statistical errors in neutron effective dose values were all less than 2%. 


z (g/cm 2 ) 

3DHZETRN (N= 1) 

3DHZETRN (N= 26) 

Geant4 

FLUKA 

PHITS 

0 

10.31 

6.98 

5.68 

5.68 

6.92 

30 

5.44 

3.50 

2.77 

2.96 

3.29 


Table 4. Neutron fluence root mean square relative differences of 3DHZETRN (N= 26) and MC codes at z = 0 
g/cm 2 (top of ICRU tissue sphere) andz = 30 g/cm 2 (bottom of ICRU sphere). Dimensionless units. 


z (g/cm 2 ) 

vs Geant4 

vs FLUKA 

vs PHITS 

MC spread 

0 

0.468 

0.308 

0.269 

0.676 

30 

0.433 

0.327 

0.339 

0.674 


Table 5. Total CPU seconds required for benchmark calculations. Total number of histories ran in Monte Carlo 
simulations are shown in pare ntheses. 


Webber SPE 

3DHZETRN (N= 26) 

59 

Geant4 

2.0x10 s (3.5x10") 

FLUKA 

5.1x10 s (3.9x10") 

PHITS 

9.6xl0 7 (1.4x10") 


Conclusions 

In this work, 3DHZETRN was extended to enable analysis of general combinatorial geometry. The 
extensions were tested by considering a cylindrical aluminum shell surrounding internal aluminum boxes and an 
ICRU tissue sphere. The internal boxes were positioned relative to the external source to provide only partial 
shielding of the ICRU sphere, thereby emphasizing 3D transport effects. Following the outline of previous studies 
[Wilson et al. 2014a-c, 2015], convergence tests were performed to determine the solution dependence on the 
number of transport directions, N. It was found that N = 26 provided nearly converged results with minimal relative 
error compared to the N= 86 solution. 

Comparisons between 3DHZETRN and the MC codes, Geant4, FLUKA, and PHITS, were made at points 
within the geometry exposed to the Webber 1956 SPE. As in the previous studies [Wilson et al. 2014a-c, 2015], it 
was found that 3DHZETRN agrees well with MC codes. Most importantly, even within the complex geometry, 
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3DHZETRN was several orders of magnitude faster than the MC simulations. Future work will provide yet another 
extension to the 3DHZETRN geometry capabilities by enabling analysis through complex and realistic geometry 
described by only a single ray-trace thickness distribution. Such geometric descriptions are typically generated in the 
vehicle design and optimization process [Wilson et al. 2003, Walker et al. 2013], where detailed computer aided 
design models are maintained and developed [Qualls et al. 2001, Simon et al. 2014]. This extension is the critical 
step in connecting and integrating 3DHZETRN into shield optimization frameworks. 
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